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Abstract 

This paper presents a semi-analytical estimate of the response of a grandstand occnpied 
by an active crowd and by a passive crowd. Filtered Ganssian white noise processes are 
nsed to approximate the loading terms representing an active crowd. Lnmped biodynamic 
models with a single degree of freedom are inclnded to reflect passive spectators occnpying 
the structure. The response is described in terms of the first two moments, employing the Ito 
formula and the state augmentation method for the stationary time domain solution. The 
quality of the approximation is compared on the basis of three examples of varying complexity 
using Monte Carlo simulation based on a synthetic generator available in the literature. For 
comparative purposes, there is also a brief review of frequency domain estimates. 

Keywords: grandstand response, random vibration, active crowd, white noise process, 
hltration 


1. Introduction 

The response of a grandstand can be resolved quite easily by linear dynamics methods, 
if we neglect all the randomness of the system. However, for a more accurate description, at 
least the most significant uncertainties need to be taken into account. The main uncertainties 
include 

• forcing terms resulting from active crowd movements, especially synchronized jumping, 

• the uncertainties of the parameters in discrete biodynamic models—randomness of 
stiffness, mass and damping matrices, 

• the size and spatial distribution of an active crowd and a passive crowd. 

Further generalizations can take into account various kinds of nonlinearities, e.g. geometrical 
and material non-linearities and non-linearities of biodynamic models (Huang and Griffin, 
2008). However, the following restrictions will be assumed for the purposes of this paper: 
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the material parameters of the structure are treated as deterministic, since their influence is 
negligible in comparison with the sources listed above and the scope of the overall response; 
the spatial distribution of the crowd is hxed; biodynamic models of the passive crowd are 
treated as deterministic. 

The results of measurements carried ont on simple structures indicate that an active 
spectator might be represented as a time dependent force, cf (Ellis and Ji, 1994). In other 
words, under certain circumstances, an active spectator does not inflnence the properties of 
the structure. Several generators of normalized artihcial load processes have been developed 
in the literature, e.g. (Sim, 2006) and (Racic and Pavic, 2010), which can be snpplemented 
by hnman body weights, e.g. (Hermanussen et ah, 2001). However, the passive part of the 
crowd is assnmed to be stationary in space and in permanent contact with the strnctnre, 
exclnding accelerations that exceed the gravity of the earth. Thus a passive spectator can be 
modeled as a biodynamic system. A survey of models of this type can be fonnd in (Sachse 
et ah, 2003). These remarks allow us to model the structure considered here in terms of 
random vibrations of linear systems. 

Monte Carlo simnlation (MC) is one of the ways used in advanced grandstand design 
procedures for reflecting all kinds of randomness. MC is a general method, bnt it has consid¬ 
erable disadvantages in its original form, e.g. rather slow convergence when estimating low 
probabilities, and high compntational demands. MC will be employed in this paper mainly 
for controlling the accuracy and the performance of semi-analytical methods introdnced in 
terms of stochastic differential eqnations. These simplihed methods can be applicable and 
usefnl in the preliminary stages of design procednres, when qnick and approximate solntions 
are snfficient. 

2. Stochastic differential equations 

Employing the Finite Element Method, and taking into acconnt the above-mentioned 
assnmptions, a mathematical model of the strnctnral system considered here can be written 
as a set of hyperbolic differential equations 

MZ{t) + CZ{t) + KZ{t) = GY{t), t>0 (1) 

where Z and Y are and M'^'-valued stochastic processes, and M, C, K and G are 
{d/2, d/2) and {d/2,d') matrices of mass, damping, stiffness and input distribution. Over¬ 
dot denotes a derivative with respect to time g{t) = d/dtg{t). To simplify subseqnent 

expressions, let us apply expectation operator E in eqnation (1) and snbstract the result 

from (1). We arrive at the set of two equations 

Milzit) + Gfizit) + K^iz{t) = G^y(f), t > 0 (2) 

MZ{t) + GZ{t) + KZ{t) = GY{t), f>0 (3) 

for the mean valne ^^^(t) = EZ(f) and centered process Z{t) = Z{t) — Hzit), f-iyi't) = 
EW(f) and Y{t) = Y{t) — f-iyit)- As will become apparent in section 3, nnder certain 
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conditions process Z (t) is approximately normal, and since differential equation (3) is linear 
with deterministic coefficients, it is reasonable to accept Gaussian approximation also for 
Y {t). This consideration leads us to stochastic differential equations and to the ltd calculus. 
Under Gaussian assumptions, the response will be completely specified by its mean 
and covariance Cz{t,s) = E[Z(t)Z{s)'^]. The required quantities can be obtained from the 
time domain or from the frequency domain. 


2.1. Solution in the time domain 

Equation (2) can be solved by direct integration or, more conveniently, by a Fourier series 
(or Fourier transform) assuming periodic mean c/ section 2.2. Let us rewrite (3) as 
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or, in a more compact form, 

X(t) = aX{t)+ bY{t), t>0 


( 4 ) 

( 5 ) 


where X is an M'^-valued state-space vector stochastic process with zero mean, and a and b 
are (d, d) and (d, d')-matrices. The solution of this differential equation is given in the form 


X(t) =0(t)X(O) + 


0{t — s)bY (s) ds 


( 6 ) 


where 0(t — s) denotes a Green function or the unit impulse response satisfying 


= a0(t - s), t > s > 0, (7) 

0(0) = I the identity and 6{t — s) = exp[a(t — s)] can be expressed as a matrix exponential, 
cf (Soong and Grigoriu, 1993). Initial conditions X(0) will be set to zero for simplicity. 
Forcing term Y[t) can also satisfy its own stochastic differential equation driven by Gaussian 
white noise W{t) = dB(t)/dt. For example, let Yiit) be a continuous-time Gaussian auto¬ 
regression scalar process of order p, denoted as AR{p), cf (Brockwell et ah, 2007). Then 
Yi(t) = efSiit) where the state vector Si(t) = [S'i 4 (t),..., 5'i,p(t)]’^ satisfies the Ito equation 


dSiit) = AiSi(t)dt + bidB(t), (8) 
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Processes of this kind are also called filtered white noise processes or colored processes, and 
they have a specific frequency content. Let us assume that Yi(t) of Y(t) = [lh(t), • • •, Yd'it)]'^ 
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are mutually independent AR{pi) processes. Then we can merge equations (4) and (8) to 
obtain one coupled system 
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( 9 ) 


where d* are column vectors with the unit in i-th position, and B{t) is an M'^'-valued Brow¬ 
nian motion. This approach is called a state augmentation method (Grigoriu, 2002). An 
extension to the case Yi{t) = J2k=i'^kit), where Yk{t) are mutually independent AR{p) pro¬ 
cesses, is carried out in an obvious manner. This methodology will be employed in section 3 
for AR{2) processes. Equation (9) can again be rewritten in compact form 

dX{t) = aX(t)dt + bdB{t), t>0, (10) 


and employing the Ito formula for semimartingales we arrive at the system of evolutionary 
equations for the response mean i-ixi't) and covariance Cx{t,s) 

aRxit)^ t > 0 , ( 11 ) 

acx{t,t) + cx{t,t)a^ + bb^, t > 0, (12) 

acx{t,s), t>s>0. (13) 

Since the driving forces dB{t) are Gaussian white noise and the coefficients are constant in 
time, the solution is an Ornstein-Uhlenbeck process with an existing stationary solution. In 
our case, stationary mean fix = 0 covariance cx(t, t) = cx(t — t) = cx = 0 which leads 
to the so-called continuous Lyapunov equation 

0 = acx + Cxa^ + bb^. (14) 


= 

Cx{t,t) = 
dcx{t,s) _ 
dt 


For details and further developments, see (Grigoriu, 2002). Since the stationary matrix Cx 
contains only response displacements and velocities, the variances of the acceleration have to 
be computed through the following formulas which are valid for weakly stationary processes 
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where denotes matrix power and denotes the stationary covariance matrix of velocities 
and accelerations. Equation (15) is evaluated employing (13), which in our special case is 
simplihed to 

exit) = Oit)cx = exp[at]cx. (16) 


2.2. Solution in the frequency domain 

Taking the Fourier transform of equation (2) leads to 

fiziu:) = if(a;)G/iy(a;) (17) 

where the FRF (Frequency Response Function) Hiu) is 

H{uj) = [-uj^M + iuC + K]-\ (18) 


i denotes a complex unit, u denotes angular frequency and gioj) denotes the Fourier transform 
of function git). Assuming Hyit) in periodic form Ay('^) = ■\/^/2 T’fc[exp{i(pfc}h(a;— 

Uk) + exp{—i(y9fc}5(a; + Ok)] where r^, (pk and Ok are amplitude, phase shift and angular 
frequency of the k-th harmonic, we obtain the mean response as a solution of the complex 
linear system of equations. 

Employing spectral decomposition of stationary random processes, the response variances 
are acquired in terms of spectral density matrices Syy nnd (‘S'yy (c<;))jj = fyiuj), where 

fyiuj) is a spectral density estimate of the centered forcing term 


/y(a;) = E / bix — u)lTiuj) dx, 


(19) 


bix) is some weight function, cf (Andel, 1976), and lyioj) denotes the corresponding peri 
odogram 
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The diagonal form of Syy suggests that we treat all input processes as independent. Knowing 
the spectral density matrix of the input vector stochastic process, we obtain the spectral 
density matrix of the output process according to Soong and Grigoriu (1993) 


SzzH = Hiu)GSyYiu)G^ Wiu) ( 21 ) 

where H^ui) denotes a Hermitian transpose to Hiu). The variance of the stationary scalar 
process Z(t) with two-sided spectral density /^(w) or with one-sided spectral density gzi^) 
is evaluated as 

poo poo 

fziuj) du= gziuj) du (22) 


4 = 


and the variance of time derivative Z(t) 
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( 23 ) 


By analogy for higher time derivatives, applying higher powers of angular frequency uj. 


5 





2.3. Transformation to modal coordinates 


It is desirable to reduce the system unknowns in equations (1), (2), (14) and (21) for 
MC, mean value, time domain and frequency domain solutions. The usual modal trans¬ 
form technique can be employed directly when a structure with proportional damping is 
loaded by active spectators only. Some problems arise when the passive part of the crowd 
is introduced. All matrices are extended according to the added degrees of freedom (dofs). 
Biodynamic models possess non-proportional damping, so complex eigenvectors should be 
used. Simultaneously, all these models have the hrst eigenfrequency close to 5 Hz (and the 
second eigenfrequency close to 8 Hz if using two-degrees-of-freedom models). This implies 
that the number of eigenvectors employed in the transform is considerably enlarged. Another 
possibility is partial modal transformation. The system matrices can be decomposed into 
four parts according to 


/ Ass Ash \ 
V Ahs Ahh J 


(24) 


where A stands either for mass, stiffness or damping matrix. Sub-matrix Ass corresponds 
to the structure, Ahh to the passive crowd, and Ash or Ahs represents the mutual in¬ 
teractions. Let us note that sub-matrix Ahh is a diagonal or a band, and Ass is sparse. 
Computing the eigenvectors corresponding to an empty structure, i.e. corresponding to Ass, 
arranged in V and writing 


V 


V 0 \ 

0 / / ’ 


(25) 


all equations can be partially transformed by V, as in the standard procedure. Then the dofs 
inherent to a passive crowd are unchanged, but the structure response is described through 
several modal coordinates. 

Concerning the Lyapunov equation (14) partially transformed into modal coordinates, 
we can further reduce the computational effort by some prior information. Let us assume 
that the system response corresponds to the equation (9) with pi = 2. Splitting all matrices 
in (14) leads to 
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where 
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-\V'^ MV]-^\V'^ KV] -\V'^MV]-'\V'^CV] 


(27) 


a 2 i = 0 and remaining sub-matrices have obvious structure. Dropped subscript at covariance 
matrix c emphasizes partial transformation to modal coordinates. Since hh^ resp. Q 22 
symmetric, in fact diagonal, the solution will be also symmetric, C 12 = c:Ji. In the case 
of AR{2) processes, sub-matrix C 22 can be computed explicitly; single AR(f2) process has 
uncorrelated state variables Si, S 2 and varS*! = 09 /( 20102 ), varS '2 = Oo/(2oi), based on 
moment equations, thus C 22 is a diagonal matrix. Introduced considerations reduce the 
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system of four equations (26) in expanded form, to the set of two coupled equations 

aiiCi2 + Ci2a^2 + «'12C22 = 0 (28) 

aiiCii + Ciiol^ + Ci2a^2 + «'i2cf2 = 0 (29) 

for unknowns Ci 2 and Cn. The set resembles Sylvester and Lyapunov equations respec¬ 
tively with reduced size. Backward transformation cx = VcV'^ gives covariance matrix for 
displacement or velocity vector, here c denotes an appropriate sub-matrix of c storing the 
modal displacements or velocities, cx then contains the nodal displacements or velocities. 


2.4- Crossings of Gaussian processes 

One of the measures of system performance is level crossing. Under some circumstances, 
in a stationary case, it can be shown that up-crossing is directly connected with the reliability 
of the system. The x-up-crossing rate of Gaussian process X{f) with non-stationary mean 
value /i(t) and stationary variance is estimated as (Soong and Grigoriu, 1993) 
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(30) 


where r'fit) is the x-up-crossing rate of level x at time t, 0(a) = 1/A/^exp —a^/2, $(«) = 
f_^(j){a) da, = varX(t), = varX(t). The total mean number of upcrossings in time 
interval [0,T] is computed according to 

ntiT)= [ (31) 

Jo 

The relations can be generalized to H-out-crossings of a d-valued stochastic process, where 
D is some set in 

Another measure of system performance from the point of view of serviceability is the 
root mean square value estimated as 


RMS 



x{tYdt 



fixity dt (j2. 


(32) 


where cx^ denotes the stationary variance of centered process X(t), and px(t) its mean value. 
Analogous formulas are valid for velocity and acceleration. 


3. Applications to the response of grandstands 

As was noted in section 1, an active spectator can be treated as a time-dependent process. 
Figure 1 shows a single realization and the spectral density of a unit process, i.e. of the 
process with Gh = 1, where Gh denotes the weight of a spectator, jumping frequency 
/ = 2.67 Hz. Spectral densities were computed from (19) with Parzen weight. The realization 
was generated according to Sim (2006). Since this function is highly periodic, we will search 
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t [s] frequency [Hz] 

(a) (b) 

Figure 1: Single time history (a) and power spectral density of 10 000 realizations (b) of 
forcing term Y{t) generated according to Sim (2006). 


the mean value in the form 

p 

fiyit) = «o + cos{k ■ 271 ft) + /3k sm{k ■ 277 ft). (33) 

k=l 

Then vector 6l of the estimated parameters do, di,..., dp, /3i,..., /3p can be found by the 
linear Least Squares Method as 

a = (34) 

where 
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ti,..., is a hne enough and equidistant partition of the time interval, f-iy = [/^y (^i), • ■ ■, 
with 

1 ^ 

^y{U) = 

' k=l 

are means over N realizations Yk{ti) in time instants tj. Generating 10 000 trajectories 
provides the coefficients given in tab. 1. Realization centered with the mean value according 
to equation (33) and the coefficients from table 1 is depicted in figure 2, together with 
the spectral density and a normalized histogram, i.e. a histogram of the process Y[t) = 
[Y{t) — /i(t)]/(jy = Y{t)/aY where ay = y/varF (t) = a/0.7627 is a stationary standard 
deviation. The centered process resembles non-Gaussian colored noise. To be more precise, 
for Y[t) we have /iy = 0 for the mean, varF = 1 for the variance, 73 y = 0.424 for the 
coefficient of skewness, and 74 y = 4.076 for the coefficient of kurtosis. For comparison, the 
standard Gaussian process has coefficients 0, 1, 0 and 3. 

Let us briefly analyze the response of a harmonic oscillator with unit mass forced by 
jumping process F (t) with / = 2.67 Hz, employing MG to justify the normality assumptions. 
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Table 1: Coefficients a for an approximation of the mean value in equation (33) for / = 2.67 
Hz. 
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Figure 2: Single centered time history of Y{t) (a), the corresponding power spectral density 
of 10 000 realizations (b), a normalized histogram with standard normal density (c) and a 
normalized histogram of a non-unit process scaled with ^/^[kN] ~ 7V(0.7709, 0.0167) (d). 
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Figure 3: Coefficient of skewness 73 and coefficient of kurtosis 74 for the normalized displace¬ 
ment and velocity of the harmonic oscillator forced by the jumping process as functions of 
eigenfrequency fi. Solid line - viscous damping ( = 0.001; dash-dot line - ( = 0.07; dashed 
line - corresponding values for the Gaussian random variable. 


The coefficients of skewness and kurtosis of the state vector X(t) = [Z(t), Z(t)]'^ as functions 
of the oscillator eigenfrequency fi for two different values of viscous damping ( are depicted 
in hgure 3. Note that for frequency range 0.5 — 7 Hz the response is approximately Gaussian. 
As was expected, worse convergence is achieved for higher damping values, cf the Rosenblatt 
theorem (Grigoriu, 1995). Normalized histograms of displacement for eigenfrequencies fi = 4 
and 12 Hz and both damping values are depicted in hgure 4. Other techniques can be 
applied for approximations of the response outside this frequency range, e.g. memoryless 
transformations of Brownian colored noise, but this lies beyond the scope of our paper. Based 
on heuristic arguments and the Gentral Limit Theorem, we can assume that the higher the 
number of active spectators, and the more complex the grandstand geometry is, the more 
Gaussian the response will be. 

A spectral density approximation of the forcing process for the frequency domain solution, 
hgure 2 (b), cannot be further simplihed. This is because the FRF of the structure has sharp 
peaks, and thus exact function values are needed. Any approximation employing indicator 
functions in the vicinities of signihcant harmonics preserving variance would be inaccurate. 

However, we can employ hltered white noise processes AR{2), which arise as a solution 
of the second order ltd equation 


C2,iYiit) + Cs^iViit) + Ci^iY,{t) ^ W(t) 


(35) 


with spectral density 


Sj(a;) 


1 

[Cl,i - 


(36) 


where (ci_j, C 2 ,i, 03 ^,) correspond to the stihness, mass and damping of a harmonic oscillator. 
This function has a sharp peak positioned at R = \/ci/c2/27r when we neglect shifts due 
to damping ehects. For a closer approximation, we assume Y{t) ~ where Yi{t) 
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(a) fi = 4 Hz, C = 0.001 (b) fi = 12 Hz, C = 0.001 



o'. . . .. . 

-5 0 5 

Data 



(c) fi = 4 Hz, C = 0.07 (d) fi = 12 Hz, C = 0.07 


Figure 4: Histograms of normalized displacement Z{t) with standard normal density based 
on 1000 MC realizations. 


are mutually independent AR{2) processes. Identification leads to a nonlinear optimization 
problem: find such Ck^i, k = 1,2,3 and i = 1,... ,n, that minimize L 2 norm || • ||l 2 of 
difference 


e(ci;, c) 


min 

c={ci,i,..., 


^ [ci,i - + (wC3,i)2 

\\e{uj,c)\\L 2 ) ca G H C M"*" compact, 

C3,n} 


(37) 

(38) 


where /y denotes a spectral density estimate of the centered force term Y{t). The problem 
can be solved by the Nonlinear Least Squares method, by Simulated Annealing etc, with 
easily estimated initial vector Cq. Optimized coefficients for n = 6 of the centered process 
Y[t) are presented in table 2 , and the corresponding spectral density and spectral distribution 
function are presented in figure 5. Note that the spectral density is two-sided, and only 
one half was integrated in the spectral distribution function, thus the variance indicated is 
0.7627/2 = 0.3814. 

The performance of the structure was briefly quantified in section 2.4, where the two 
expressions (30) and (32) depended on the response variance. Thus an alternative approach 
is to optimize the response variance directly across some eigenfrequency range of a harmonic 
oscillator. Such coefficients are summarized in table 2, with indices var, spectral density and 
distribution function are depicted in hgure 5, frequency range 0.5 — 10 Hz. 

Let us briefly note the situation when the forcing process Y[t) is not a unit process, 
i.e. Gh 7 ^ L Since we are limiting our considerations to Gaussian approximation, only the 
first two moments of Gh will apply. The deterministic weight corresponds to a singular 
case yarGn = 0. Then the forcing term has the form Ycit) = GnYlt), mean response 
l^z,G{t) = ^[GH]fJ'z(t) and stationary response variance cr|g = E[G|^](t|, where Hzif) and 
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Table 2: Coefficients Ck^i, k = 1,2,3 and i = 1,..., 6 of the six independent AR{2) members 
used for approximation of the centered forcing term Y{t) in frequency range 0.5 — 10 Hz. 


i 


C2,i 


\/ci,i/c2,i/27r 

1 

90.8657 

0.3227 

0.0148 

2.67 

2 

35.9464 

0.1276 

0.1167 

2.67 

3 

283.6701 

0.2520 

0.0118 

5.34 

4 

74.1544 

0.0661 

0.0737 

5.33 

5 

913.9890 

1.1120 

21.5186 

4.56 

6 

228.5270 

0.0907 

0.1576 

7.99 

Ivar 

91.0909 

0.3237 

0.0076 

2.67 

2 

40.3066 

0.1430 

0.0804 

2.67 

3var 

281.0462 

0.2490 

0.0209 

5.35 

4:var 

83.7399 

0.0746 

0.0573 

5.33 

5var 

914.0015 

0.9376 

21.6921 

4.97 

6 

'-^var 

228.7642 

0.0908 

0.1552 

7.99 




(a) spectral density 


(b) spectral distribution function 


Figure 5: Spectral density and spectral distribution function of centered forcing Y(t) and of 
its approximation where Yi{t) are independent AR{2) processes with the coeffi¬ 

cients in table 2 based on spectral and variance optimization. 
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Data Data Data 


Figure 6: Tested structures, geometries with labeled points of interest and normed histograms 
of the displacement, structures occupied by an active crowd only. 

Table 3: First eigenvalues corresponding to vertical bending modes [Hz], cf figure 6. 


Structure 

fi 

f2 

fa 

f4 

fs 

fe 

f? 

fs 

beam 

7.5 

22.5 

40.3 

— 

— 

— 

— 

— 

cantilever 

5.4 

7.0 

8.2 

25.1 

28.4 

— 

— 

— 

grandstand 

2.5 

2.6 

3.8 

4.5 

4.8 

4.9 

5.4 

6.0 


ct| are the response mean and the variance of the structure loaded by the unit forcing term 
Y{t). However we should be aware that even when process Y{t) was a Gaussian, process 
Ycit) as a product of a random variable with a stochastic process, is not Gaussian. To 
quantify the influence of such scaling, compare the histograms in hgure 2 (c) and 2 (d), 
where Gh has normal distribution with mean value 0.7709 and variance 0.0167. 

4. Numerical examples and comparison 

In this section, the quality of the approximation in the time or frequency domain will 
be compared with MG simulation. A total of four mechanical systems will be tested: a 
harmonic oscillator, a simply supported beam, a simple cantilever grandstand, and a realistic 
grandstand with 1, 4, 72 and 630 positions for spectators, respectively. Geometries with 
centered normed response histograms for an active crowd only, based on MG simulation, 
are depicted in hgure 6, and several lowest eigenfrequencies corresponding to the vertical 
bending modes are presented in table 3. Note also that assumptions on convergence to 
normal distribution are approximately fulhlled. All examples are artihcial, not realistic, so 
the measured responses would be unacceptable. 
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Figure 7: Response mean displacement in comparison with an approximation based on the 
hrst four harmonics for a harmonic oscillator, /i = 5 Hz ^ = 0.07. 

4.1. Harmonic oscillator 

Let us assume a harmonic oscillator with unit mass, two values of viscous damping ( 
and variable stiffness. The mean value response is presented in hgure 7. The response was 
acquired by direct integration of (2), starting at t = 0.5 s. The approximation utilizes 
equation (33) and the coefficients in table 1. A comparison of the total mean up-crossings 
n+(T) in the time interval [0,T], T = 160 s, as functions of oscillator eigenfrequency fi for 
two values of viscous damping ( = 0.001 and ( = 0.07 and for two hxed levels x = 0.002 
and X = 0.005 m are depicted in hgure 8, employing formulas (30) and (31). The size of 
time interval T is based on heuristic considerations about the average length of the musical 
compositions. The time domain solution is based on the sum of six independent AR{2) 
processes with the coefficients in table 2 for var optimization. The stationary response 
variances are computed according to formulas (28), (29) and (15). The frequency domain 
approximation employs equations (21), (22) and (23). The spectral density estimate /y (w) of 
the centered input process has 307 values over the frequency range 0 — 10 Hz, using variable 
division. Figure 8 shows that the results are roughly in agreement with MC in the frequency 
range 0.2 — 10 Hz. Note that the total number of mean zero-up-crossings for a deterministic 
periodic function with frequency 2.67 Hz is 160 ■ 2.67 ~ 427, cf hgure 8 (c) and (d), where 
distinct plateaux are found. The results for MC are based on 1000 realizations 160 s in 
length. 

4-2. A simply supported beam 

The next example is a simply supported beam with Rayleigh damping = 0.05 and 
C ,2 = 0.08 for the hrst two vertical modes and total mass 1700 kg. Poor approximations are 
anticipated, since the structure is quite stih with high eigenfrequencies, see table 3 and the 
results for the harmonic oscillator. Two cases are studied: a structure occupied by an active 
crowd only; and a structure occupied by a mixed crowd. In the second case, the two left 
hand side positions are loaded by forces and the two right hand side positions are occupied 
by passive spectators. Deterministic biodynamic models according to Coermann are used. 
These are simple-degree-of-freedom oscillators with mass 86.2 kg, stihness 85.25 kN/m and 
viscous damping 1.72 kNs/m, eigenfrequency 5 Hz, cf (Sachse et ah, 2003). Results are 
presented only for the labeled point in hgure 6. The total mean up-crossings for the hrst case 
n+(160) as a function of level x are depicted in hgure 9 (a). The RMS values for acceleration 
are 1.947 m/s^ for the MC solution, 2.135 m/s^ for the frequency domain solution, and 
1.914 m/s^ for the time domain solution. All input processes are treated as independent, so 
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(a) C = 0.001 



(b) C = 0.001 




(c) C = 0.07 


(d) C = 0.07 




(e) C = 0.001 


(f) C = 0.07 


Figure 8: Total mean up-crossings n+(T) and acceleration RMS values as functions of fi for 
a harmonic oscillator, distinct levels x and viscous damping T = 160 s. 
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(a) active crowd only (b) 2 active, 2 passive spectators 
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Figure 9: Total mean up-crossings u+(T) of a simply supported beam as functions of x, 
T = 160 s for an active crowd (a) and for a mixed crowd (b), single realization for an active 
crowd (c) and mean response for an active crowd (d). 


matrix Syy equation (21) has nonzero only diagonal entries. The mean value response 
is depicted in hgure 9 (d) with a single realization (c). The total up-crossings for a mixed 
crowd are depicted in hgure 9 (b), and the RMS values are 1.048 m/s^ for the MC solution, 
1.094 m/s^ for the frequency domain solution, and 0.990 m/s^ for the time domain solution. 
Let us also recall our assumption of hxed spatial distribution of the crowd, mass coefficient 
7 = mu/'ms = 0 . 1 , where mn denotes the total mass of passive spectators, and ms denotes 
the total mass of the structure. Results for MC based on 2000 realizations. The approximate 
shape of the total up-crossings, especially for a mixed crowd, differ from the MC simulation, 
because of the non-Gaussian response due to high structure eigenfrequencies. 

4 . 3 . Cantilever grandstand 

This system has total mass 18.2 t, geometry according to hgure 6 , and is loaded with 72 
active spectators in the hrst case, Rayleigh damping with (^i = 0.05 and (2 = 0.08 is used for 
the hrst two vertical modes. The total up-crossings of the response displacement are depicted 
in hgure 10 (a), and single realization and the mean response are depicted in sub-hgures (c) 
and (d), RMS accelerations 5.583 m/s^ for the MC solution, 5.682 m/s^ for the frequency 
domain solution, and 5.616 m/s^ for the time domain solution. For 36 spectators chosen to 
be passive according to Coermann with uniformly random but hxed positions, the resulting 
up-crossings are depicted in sub-hgure (b), mass coefficient 7 = 0.17. The acceleration RMS 
values in this case appear to be 1.327 m/s^ for the MC solution, 1.347 m/s^ for the frequency 
domain solution, and 1.309 m/s^ for the time domain solution. The results for MC are again 
based on 2000 realizations 160 s in length. 
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(a) active crowd only (b) 36 active, 36 passive spectators 




t[s] 


(d) response mean 


Figure 10: Total mean up-crossings of a cantilever grandstand, active crowd (a) and mixed 
crowd (b), single realization (c) and mean response (d) for an active crowd. 


4-4- Realistic grandstand 

In this concluding example, let us briefly examine the results acquired for a complex 
structure. The geometry is sketched in hgure 6 with the point of interest labeled, total mass 
148.6 t, Rayleigh damping with (^i = 0.01 for the first and (2 = 0.02 for the sixth vertical 
mode used here. The results for the structure loaded by an active crowd only are depicted 
in hgure 11 (a) (c) (d). The acceleration RMS values appear to be 5.546 m/s^ for the MC 
solution, 5.565 m/s^ for the frequency domain solution, and 5.597 m/s^ for the time domain 
solution. In the second case, 315 positions are occupied by passive spectators according 
to Coermann, mass coefficient 7 = 0.18. The results for this case are presented in hgure 
11 (b), RMS accelerations 1.905 m/s^ for the MC solution, 1.904 m/s^ for the frequency 
domain solution, and 1.919 m/s^ for the time domain solution. MC simulation based on 
2000 realizations. 

4-5. Comparison and performance 

The time consumption for the diherent solution techniques is summarized in table 4, 
where the demands of MC simulation are presented for 50 realizations only. This value 
is based on the convergence tests of total up-crossings presented in hgures 12 (a) and (b) 
for a cantilever grandstand. Obviously these tests are highly sensitive to level x, and the 
value used can be considered as a lower bound. The size of the time integration step is 
chosen to be h = 0.01 s, Newmark integration scheme used. The number of dofs of each 
system is also mentioned, together with the size of the Lyapunov equation (14) ul for 
a mixed crowd. Obviously ul = 2ndof + 2 np -|- 12na where Udof denotes the number of 
dofs of the empty structure, Up is the number of dofs of the passive spectators, and is 
the number of active spectators. In the case of a cantilever grandstand we therefore have 
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(a) active crowd only 


(b) 315 active, 315 passive spectators 
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Figure 11: Total mean up-crossings of a realistic grandstand, active crowd (a), mixed crowd 
(b), single realization (c) and mean response (d) for an active crowd. 



Figure 12: MC convergence tests of n+(160) for a cantilever grandstand, cf figure 10 (b). 

2 ■ 504 -|- 2 ■ 36 -|- 12 ■ 36 = 1512, since passive spectators are modeled as systems with a single 
degree of freedom. For the purposes of comparison, the time consumption for the partial 
modal transform are also summarized together with the number of eigenvectors used Ueig 
and according highest frequency feig. All simulations were performed on core i7 with a 16 GB 
RAM computer, Matlab® parallel implementation. The results in the table suggest that 
when 50 realizations are employed, only structures of moderate size can be effectively solved 
by semi-analytical methods. 

5. Conclusions 

This paper has presented a study of the vibration of grandstands loaded by an active 
crowd, using Gaussian approximation of the response. The main results can be summarized 
as follows: 

1 . A mathematical description of the response of a mechanical system employing spectral 
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Table 4: Comparison of computational demands. 


Method/Struct. 

Harm. osc. 

Beam 

Cant, grand. 

Real, grand. 

full system 

'^dof 

1 

29 

504 

4 068 

Ul 

14 

86 

1512 

12 546 

MC 50 

0.601 s 

20.078 s 

380.451 s 

5 484 s 

Freq. domain 

0.106 s 

2.672 s 

74.639 s 

4 178 s 

Time domain 

0.112 s 

2.392 s 

38.915 s 

6 507 s 

partial modal transform, cf equations (27), (28) and (29) 

Tleig 

~ 

7 

10 

50 

f . 

ieig 

~ 

37 Hz 

30 Hz 

23 Hz 

ul 

~ 

42 

524 

4 510 

MC 50 

- 

13.240 s 

36.214 s 

307 s 

Freq. domain 

- 

1.764 s 

4.452 s 

734 s 

Time domain 

- 

1.588 s 

2.582 s 

67 s 


and time domain solutions for weakly stationary Gaussian excitations has been recalled. 
Partial modal transformation due to a passive crowd has been briefly discussed. 

2 . A motivating example of a harmonic oscillator has shown that the normalized dis¬ 
placement and velocity have approximately normal distribution under the conditions 
on eigenfrequencies and damping. 

3. Taking this fact into account, the mean value of the force has been approximated as 
a truncated Fourier series, and the spectral density of the centered process has been 
estimated employing the Parzen window for the frequency domain solution. For the 
time domain solution, we have employed a linear combination of independent auto¬ 
regression processes of the second order with coefficients optimized to achieve the least 
error in the response variance of a harmonic oscillator. 

4. Three different examples of varying complexity have shown the quality of the response 
approximation in terms of total displacement up-crossings and acceleration RMS in 
comparison with Monte Carlo simulation. Limitations following from a simple oscillator 
have been confirmed on multi-degree-of-freedom systems. 

5. The computational demands have been measured and summarized in terms of the time 
needed for solution, and the applicability of the techniques has been proved. 

Finally, let us note that these methods are approximations and can be further refined. In 
particular, the solution can be reformulated for the non-Gaussian processes for which higher 
moments can be derived. These apply in improving the mean up-crossing rate estimates for 
stiff structures. 
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